// /opt/openmpi/bin/mpic++ -lm -o m mpiden.cpp
#include <mpi.h>
//    #undef SEEK_SET
//    #undef SEEK_END
//    #undef SEEK_CUR

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <queue>
//#include <procfs.h>
#include "/usr/include/sys/procfs.h"

using namespace std;
#include "general.h"
#include "mpiden_objects.h"
#include "mpiden_serialvars.h"
#include "mpiden_operations.h"
#include "mpiden_functions.h"

// masa por particula de la sim de 500Mpc es  7.1973e+10 masa solares
// masa por particula de la sim de 60Mpc  es  1.056e+8 masa solares

#define WRITERS 4      // Numero de files que se escribiran en el output
#define INTBOX   500   // 60
#define SNAPNUM  37    // 24
#define FILES    8     // numero de files por snapshot 
#define NPAQMAX 12000  // numero maximo de particulas por paquete de comunicacion

// Structuras a comunicar
MPI_Datatype  Dosint;
MPI_Datatype  Partmpi; 

int  totalsizeR,totalsizeL,maxsizeR,maxsizeL;

struct equival *arista   ;            // conjunto de equivalencia en grupos con el nodo de la izquierda

int main(int argc, char **argv)
{
  int i,j,ibox;	
  int iy, iz ,sendproc, recvproc;
  double numdens;
  double r0;
  FILE *pfout;
  long int kkk;
  char fileout[200];
  MPI_Datatype  type[2]= {MPI_INT, MPI_INT};
  int blocklen[2]      = { 1, 1};
  MPI_Aint disp[2];
  
  MPI_Init (&argc, &argv);
  MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank (MPI_COMM_WORLD, &proc);

  // Crea la estructura Dosint
  interDer  = new struct face [2] ;
  disp[0] = (MPI_Aint)&interDer[0].Id - (MPI_Aint)&interDer[0];
  disp[1] = (MPI_Aint)&interDer[0].Grp - (MPI_Aint)&interDer[0];
  MPI_Type_create_struct(2, blocklen, disp, type, &Dosint);
  MPI_Type_commit(&Dosint);
  delete[] interDer;

  read_header();

  numdens=(double)nparticles_total/pow(boxsize,(double)3.); 
  r0 = (double)0.17*pow(numdens,(double)(-1./3.));
  overhead = r0/2.;
  cout << overhead << endl;
 // MPI_Barrier(MPI_COMM_WORLD);
 // FIN;

    
  nobj   = read_gadget_overlap();

//  cout << "despues de leer" << endl;  fflush(stdout);
  xsize = (boxsize/(double)(numprocs-1) + (double)2.*overhead);

  np=(int)floor(boxsize/overhead); if(np>512)np=512 ; //
  nx=(int)floor(xsize  /overhead); if(nx>512)nx=512 ;
  
  facbinx= 1./overhead;//(double)nx/xsize   ;
  facbinp= (double)np/(double)boxsize ;

//  printf("Box %f Proc %d / %d\n",boxsize,proc,numprocs);fflush(stdout);
  

  if(proc<numprocs-1)
  {


  change_coord();
  
  
 /* pfout=fopen("/dev/urandom","r");
  fread(&kkk,sizeof(long int),1,pfout);
  fclose(pfout);
  srand48(kkk);
   
  sprintf(fileout,"Particles.%d.dat",proc);
  pfout=fopen(fileout,"w");
  if(pfout==NULL) {moco("Open failed\n");FIN;}

  for(i=0;i<particle.size();i++)
  {
	  if(drand48()<0.1)
	  {
             fprintf(pfout,"%f %f %f %f %f %f %d\n",
	           particle[i].Pos[0],particle[i].Pos[1],particle[i].Pos[2],
                   particle[i].Vel[0],particle[i].Vel[1],particle[i].Vel[2],
	           particle[i].Id);
	  }
  }*/
  

  RED; 
  cout << " proceso " << proc ;  
  CYAN;
  cout << " r0 "<< r0 ;
  cout << " np "<< np ;
  cout << " nobj "<< nobj ;
  cout << " nx "<< nx << endl;  DEFA;
  
  try
  {
	  lirst      = new int [nx*np*np];
	  ll         = new int [particle.size()];
  }
  catch(bad_alloc) {cout << "la cago en el lirst y ll"<< flush << endl;}

  //fclose(pfout);
  }

  grid();
//  cout << totalsizeR << "  |||  " << totalsizeL << endl;
//   FIN; 
//

  MPI_Barrier(MPI_COMM_WORLD);
  
  idenBFS((float)r0);
  if(proc==numprocs-1)ngrupos=0;
  MPI_Allreduce(&ngrupos,&maxngrp,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD); 

  //////
  MPI_Barrier(MPI_COMM_WORLD);
  cout << "en 1" << proc << endl;

  interDer  = new struct face [maxovh*2] ;
  interIzq  = new struct face [maxovh*2] ;
  externIzq = new struct face [maxovh*2] ;
  for(i=0;i<maxovh;i++)
  {
	  interDer[i] .Id=-1 ;
	  interIzq[i] .Id=-1 ;
	  externIzq[i].Id=-1;
	  interDer[i] .Grp=-1;
	  interIzq[i] .Grp=-1;
	  externIzq[i].Grp=-1;
  }

  interface(interDer,nx-1);// check
  interface(interIzq,0);// check

  MPI_Barrier(MPI_COMM_WORLD);
  cout << "en 2 " << proc << endl;

  if(proc<numprocs-1)
  {
	  sendproc = proc; recvproc = proc;
	  if(proc < numprocs - 2)
		  sendproc = proc + 1 ;
	  else
		  sendproc = 0;
	  if(proc > 0)
		  recvproc = proc - 1 ;
	  else
		  recvproc = numprocs - 2 ;

	  cout<< proc << " | "
		  << sendproc << " | "
		  << recvproc << endl;


	  if(proc%2==0)
	  {
	      	  MPI_Send(interDer, maxovh,Dosint,sendproc,34,MPI_COMM_WORLD);
	      	  MPI_Recv(externIzq,maxovh,Dosint,recvproc,53,MPI_COMM_WORLD,&status);
	  }
	  else
	  {
	      	  MPI_Recv(externIzq,maxovh,Dosint,recvproc,34,MPI_COMM_WORLD,&status);
	      	  MPI_Send(interDer, maxovh,Dosint,sendproc,53,MPI_COMM_WORLD);
	  }
  } 
  else 
  {  
	  delete[] interIzq ;  
	  delete[] externIzq;
  }
  
  delete[] interDer;

  MPI_Barrier(MPI_COMM_WORLD);
  cout << "en crrearista " << proc << endl;
  crearista();// check
  cout << "fuera de crrearista " << proc << endl;
  
  if(proc<numprocs-1)
  {
  delete[] interIzq;  delete[] externIzq;
  }
  ///////

  MPI_Barrier(MPI_COMM_WORLD);
  cout << "en creagrafo " << proc << endl;
  creagrafo();// check
  MPI_Barrier(MPI_COMM_WORLD);
  cout << "fuera de creagrafo " << proc << endl;

  //numero   = new int [maxngrp*numprocs-1,0];
  //totalnum = new int [maxngrp*numprocs-1,0];
  
  MPI_Barrier(MPI_COMM_WORLD);
  cout << "antes de names " << proc << endl;
  groups_namesBFS();// check
  cout << "despues de names " << proc << flush << endl;
  MPI_Barrier(MPI_COMM_WORLD);
  MPI_Finalize ();
  exit(0);
    
  MPI_Allreduce(&numero,&totalnum,(maxngrp*numprocs-1),MPI_INT,MPI_SUM,MPI_COMM_WORLD);
  //delete [] numero;
  
  write_data(WRITERS);// check
  
  MPI_Finalize ();
  exit(0);
}

int write_data(int nfiles)
{
	int i,h,j,k;
	int recvproc, filegrp;
	int comm_map[numprocs];
	int send_map[numprocs];
	struct ptcl *commbuff;
	vector <struct ptcl> halo;
        ofstream Salida ;
	char *filename=NULL;

	filegrp =(int) ceilf((float)universalname/(float)nfiles);

        sprintf(filename,"Salida.%d.dat",proc);
	if(proc<nfiles) Salida.open(filename);

	for(recvproc=0;recvproc<nfiles;recvproc++)
	{
	     for(i=recvproc*filegrp; i<(recvproc+1)*filegrp -1 && i<universalname;i++)
	     {
		     for(j=0;j<numprocs;j++){ comm_map[j]=0; send_map[j]=0;}
		     comm_map[proc]=numero[i];
                     MPI_Allreduce(&comm_map,&send_map,numprocs,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
		     if(proc==recvproc)
		     {
			     for(h=0;h<numprocs;h++)
			     {
				     if(send_map[h]>0)
				     {
					     commbuff= new struct ptcl [send_map[h]];
                                             MPI_Recv(&commbuff,send_map[h],Partmpi,h,34,MPI_COMM_WORLD,&status);
					     for(k=0;k<send_map[h];k++)
						     halo.push_back(commbuff[h]);
					     delete[] commbuff;
				     }
				     
			     }
			     if(halo.size()!=totalnum[i]){moco("la cagaste al escribir\n"); FIN;}

                             // Escribiendo la salida
			     Salida << halo ;
		     }
		     else
		     {
			     if(comm_map[proc]>0)
			     {
				     commbuff= new struct ptcl [comm_map[proc]];
				     for(h=0;h < grupo[i].members.size();h++)
					     commbuff[h]=*(grupo[i].members[h]);
			     	     MPI_Send(&commbuff,comm_map[proc],Partmpi,recvproc,34,MPI_COMM_WORLD);
				     delete[] commbuff;
			     }
		     }

	     }
	}
  
	if(proc<nfiles) Salida.close();
	return 0;
}

// modificada para que sobre un micro
int groups_namesBFS(void)
{
	int i,j,*rename,sendproc,h;
	int *listR,*listL;
	int *HeadR,*HeadL;
	long int kk;
	float percent;
	struct point u,v,der,izq;
	vector <point>::iterator ptr ;
	cola Q;

	MPI_Barrier(MPI_COMM_WORLD);
	cout << maxngrp << " estamos? 1 " << proc << endl;
	try{rename = new int [maxngrp];}catch(bad_alloc){moco("rororororor \n");FIN;}
	rename[700000]=1;
	cout << proc << " " << rename[700000] << " ---------" << maxngrp << endl;
	if(proc==numprocs-1)
	{

	   try {HeadR = new int [grafo_vertices.size()];}catch(bad_alloc){moco("kk1\n");FIN;}
           try {HeadL = new int [grafo_vertices.size()];}catch(bad_alloc){moco("kk2\n");FIN;}
           try {listR = new int [grafo_vertices.size()];}catch(bad_alloc){moco("kk3\n");FIN;}
           try {listL = new int [grafo_vertices.size()];}catch(bad_alloc){moco("kk4\n");FIN;}
	   for(i=0;i<grafo_vertices.size();i++) { HeadR[i]=-1; HeadL[i]=-1;}
	   for(i=0;i<grafo_vertices.size();i++) 
	   { 
		   if((*(grafo_arista[i].r)).name>grafo_vertices.size()-1){moco("hhhhhh "); FIN;}
		   if((*(grafo_arista[i].l)).name>grafo_vertices.size()-1){moco("llllll "); FIN;}
		   HeadR[(*(grafo_arista[i].r)).name]=i; 
	     	   HeadL[(*(grafo_arista[i].l)).name]=i;
	   }
		MPI_Barrier(MPI_COMM_WORLD); cout << "llegaroto !!" <<endl;
	FIN;
	   for(i=0;i<grafo_vertices.size();i++) 
	   { 
		   listR[HeadR[(*(grafo_arista[i].r)).name]]=i; 
		   HeadR[(*(grafo_arista[i].r)).name]=i; 
	     	   listL[HeadL[(*(grafo_arista[i].l)).name]]=i;
	     	   HeadL[(*(grafo_arista[i].l)).name]=i;
	   }
		   
	   cout << "conectando grupos" << flush << endl;
	   for(i=0;i<grafo_vertices.size();i++)
	   {   
	       percent= ((float) i/(float)grafo_vertices.size())*100.;
	       if(fmodf(percent,1.)==0.)
	    	       cout << "conectando "<< percent  << " % : proc: " << proc << endl;
	       if(fmodf((float)i,1000.)==0.)
		       cout << i << " / " << grafo_vertices.size() << '\r' ;
	       
               if(grafo_vertices[i].visitado) continue ;
               grafo_vertices[i].visitado  = true ;
               grafo_vertices[i].distancia = 0    ;
               //Encolar(Q, s);
	       Q.push(grafo_vertices[i]);
               do
               {
                  // extraemos el nodo u de la cola Q y exploramos 
	          // todos sus nodos adyacentes
                  u     = Q.front();
   	          Q.pop();
	       
                  //for(j=0;j<grafo_arista.size();j++)
		  /*
                  for(j=0;j<10;j++)
	          {
		       if( u.name == (*(grafo_arista[j].r)).name)
		       {
			     if(!(*(grafo_arista[j].l)).visitado)
			     {
                                  (*(grafo_arista[j].l)).visitado    = true            ;
                                  //(*(grafo_arista[j].l)).distancia   = u.distancia + 1 ;
		                  //(*(grafo_arista[j].l)).padre = &grafo_vertices[(*(grafo_arista[j].r)).name];
		                  (*(grafo_arista[j].l)).padre = grafo_vertices.begin() + (*(grafo_arista[j].r)).name ;
		                  //(*(grafo_arista[j].l)).padre = find_if( grafo_vertices.begin(), grafo_vertices.end(), 
				  //				look_vertice_name(*(grafo_arista[j].r)) );
	                          Q.push(*(grafo_arista[j].l))     ;  //Encolar v en Q
			     }
		       }
		       if( u.name==(*(grafo_arista[j].l)).name)
		       {
			     if(!(*(grafo_arista[j].r)).visitado)
			     {
                                  (*(grafo_arista[j].r)).visitado    = true            ;
                                  //(*(grafo_arista[j].r)).distancia   = u.distancia + 1 ;
		                  //(*(grafo_arista[j].r)).padre       = &grafo_vertices[(*(grafo_arista[j].l)).name];
		                  (*(grafo_arista[j].r)).padre = grafo_vertices.begin() + (*(grafo_arista[j].l)).name ;
		                  //(*(grafo_arista[j].r)).padre = find_if( grafo_vertices.begin(), grafo_vertices.end(), 
				  //			look_vertice_name(*(grafo_arista[j].l)) );
	                          Q.push(*(grafo_arista[j].r));//Encolar(Q, v);
			     }
		       }
	          }

		  */
	
	///	  
		  h=HeadR[u.name];
		  if(h!=-1)
		  {
			  do
			  {
				  h=listR[h];
				  if(!(*(grafo_arista[h].l)).visitado)
				  {
				       	  (*(grafo_arista[h].l)).visitado    = true            ;
				     	  (*(grafo_arista[h].l)).padre = grafo_vertices.begin() + (*(grafo_arista[j].r)).name ;
				     	  Q.push(*(grafo_arista[j].l))     ;  //Encolar v en Q
			       	  }
				  if(h==HeadR[u.name])break;
			  }while(1);
		  }
		  
		  h=HeadL[u.name];
		  if(h!=-1)
		  {
			  do
			  {
				  h=listL[h];
				  if(!(*(grafo_arista[h].r)).visitado)
				  {
				       	  (*(grafo_arista[h].r)).visitado    = true            ;
				     	  (*(grafo_arista[h].r)).padre = grafo_vertices.begin() + (*(grafo_arista[j].l)).name ;
				     	  Q.push(*(grafo_arista[j].r))     ;  //Encolar v en Q
			       	  }
				  if(h==HeadL[u.name])break;
			  }while(1);
		  }
		  
		  
               }while(!Q.empty());
	       //agotadas todas las adyacencias salimos
	   }
	   
           delete [] HeadR;
           delete [] HeadL; 
           delete [] listR; 
           delete [] listL; 

	   cout << "renombrando" << endl; 

	   universalname=0;
	   for(i=0;i<grafo_vertices.size();i++)
	   {
		u=grafo_vertices[i];
	        ptr=u.padre;
	
	        // es padre ?
		if( ptr ==(vector <point>::iterator )NULL)
		{
			grafo_vertices[i].name=universalname;
			universalname++;
		}
	   }

           MPI_Bcast(&universalname,1,MPI_INT,numprocs-1,MPI_COMM_WORLD);

	   for(i=0;i<grafo_vertices.size();i++)
   	   {
		u=grafo_vertices[i];
	        ptr=u.padre;
		
		// es hijo de
		if( ptr !=(vector <point>::iterator )NULL)
		{
			
			do
		 	{
			     v=*ptr; ptr=v.padre;
                             if(ptr==(vector <point>::iterator )NULL)
				     break;
  		        }while(1);
	   	        grafo_vertices[i].name = v.name ;
	   	}
 	   }
	   
	   sendproc=0;h=0;
	   for(i=maxngrp;i<grafo_vertices.size();i++)
   	   {
		   rename[h++]=grafo_vertices[i].name;
		   if(h==maxngrp)
		   {
//			MPI_Send(rename,maxngrp,MPI_INT,sendproc++,66,MPI_COMM_WORLD);
			if(sendproc>numprocs-1)moco("la cagaste en groups_namesBFS\n");
			h=0;
		   }
	   }

           for(i=0;i<maxngrp;i++)
		   rename[i]=grafo_vertices[i].name;
	}
	else
	{
		MPI_Barrier(MPI_COMM_WORLD); cout << "llegamos caca" <<endl;
//	FIN;
//		MPI_Recv(rename,maxngrp,MPI_INT,numprocs-1,66,MPI_C/OMM_WORLD,&status);
	}

	MPI_Barrier(MPI_COMM_WORLD);
	cout << "estamos? 2 " << proc << endl;

	if(proc<numprocs -1)
	{
	for(i=0;i<particle.size();i++)
	{
	        particle[i].Idgrupo = particle[i].Idgrupo>=0 ? maxngrp*proc + particle[i].Idgrupo : -1 ;//BUG2
		particle[i].Idgrupo = rename[particle[i].Idgrupo];
                if(particle[i].Idgrupo>=0)numero[particle[i].Idgrupo]++;
	}
        }
	
	return 0;
};

// modificada para que sobre un micro
int creagrafo(void)
{
  int i,j;
  struct link  tmplink;
  struct point tmppoint,look;

 if(proc==numprocs-1)
 {
	 cout << "Hello! soy el root!" << endl;
	 for(i=0;i<maxngrp*(numprocs-1);i++)
       	 {
		 tmppoint.name = i                                 ;
		 tmppoint.visitado = false                         ;
		 tmppoint.distancia = 0                            ;
	      	 tmppoint.padre = (vector <point>::iterator )NULL  ;
	       	 grafo_vertices.push_back(tmppoint)                ;
	 }
	
	 for(i=0;i<numprocs-1; i++)
	 {
             cout << maxequival << "asdasd  "<<endl;
             MPI_Recv(arista,maxequival,Dosint,i,40,MPI_COMM_WORLD,&status);
             cout << "hola!!!!!!!!!!!  "<< i <<endl;

	     for(j=0;j< maxequival;j++)
	     {
	         if(arista[j].local!=-1 && arista[j].externa!=-1)
		 {
		     //look.name=arista[j].local;
       		     tmplink.r = grafo_vertices.begin() + arista[j].local;
       		     tmplink.l = grafo_vertices.begin() + arista[j].externa;
		     
		     grafo_arista.push_back(tmplink) ;
     	             if(arista[j].local>maxngrp*numprocs-1 || arista[j].externa>maxngrp*numprocs-1){cout <<"mal 100"; FIN;}
	         }
	     }
	     
	 }
  } 
  else
  {
	  cout << maxequival << "XXXXX  "<<endl;
     	  MPI_Send(arista,maxequival,Dosint,numprocs-1,40,MPI_COMM_WORLD);
  }
  
  return 0;
}

// modificada para que sobre un micro
int crearista(void)
{
  int i,j;
  vector <struct equival> palito;
  vector <equival>::iterator older;
  struct equival ver;
  int idcheck, sisa;

  // cargando las aristas locales
  if(proc<numprocs-1)
  {
	  for(i=0;i<maxovh;i++)
	  {
	     	  ver.local=interIzq[i].Grp;
	     	  if(ver.local==-1) continue ;
	     	  idcheck=interIzq[i].Id;
	     	  ver.externa=-2;
	     	  for(j=0;j<maxovh;j++)
		  	  if(externIzq[j].Id == idcheck)
				  ver.externa=externIzq[j].Grp;
		  /*
	     	  if(ver.externa==-2) 
		  { 
			  moco("la cagaste en crearista ") ; 
			  cout << proc << " ovh: " 
			  << maxovh << " numeros: " 
			  << ver.local <<" | " 
			  << ver.externa<< endl; 
		  }*/
		  if(ver.externa==-2) continue ;
		  if(ver.externa==-1) continue ;
		  older = find_if( palito.begin(), palito.end(), 		
			     	  comp_arista(ver) );
	     	  if (older != palito.end()) continue ;
	       	  palito.push_back(ver);
	  }
 	  sisa = palito.size();
  } 
  else 
  {
       	  sisa=0;
  }

  MPI_Allreduce(&sisa,&maxequival,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
  arista = new struct equival [maxequival] ;


  if(proc<numprocs-1)
  {
  for(i=0;i<maxequival;i++)
  {
      if(i<palito.size())
      {
          arista[i].local   = palito[i].local;
          arista[i].externa = palito[i].externa;
      }
      else
      {
          arista[i].local   = -1;
          arista[i].externa = -1; //BUG pa arreglar
      }
  }
  palito.clear();
  }
  return 0;
};

int interface(struct face *franja, int ix)
{
  int i, h,ibox;
  int iy, iz;
  char fileout[100];
  
  
  if(proc==numprocs-1)return 1;

//  ofstream Salida ;
//  sprintf(fileout,"interface.%d.proc_%d.dat",ix,proc);
//  Salida.open(fileout);
// cout << proc << " interface <<"<< endl ;
  h=0;
  for( iy = 0 ; iy < np ; iy++)
  { 
      for( iz = 0 ; iz < np ; iz++)
      {
        ibox=(ix * np + iy) * np + iz ;
        i=lirst[ibox];
        if(i==-1)continue;
        do
        {
          i=ll[i] ;
	  
	  if(((double)particle[i].Pos[0]         < 2.*overhead && ix == 0   ) ||
	     (xsize - (double)particle[i].Pos[0] < 2.*overhead && ix == nx-1)    )
	  {
//	     Salida << particle[i].Id << " ";
//	     Salida << particle[i].Pos[0] << " ";
//	     Salida << particle[i].Pos[1] << " ";
//	     Salida << particle[i].Pos[2] << " ";
	     franja[h].Id  = particle[i].Id ;
	     franja[h].Grp = particle[i].Idgrupo>=0 ? maxngrp*proc + particle[i].Idgrupo : -1 ;
//	     Salida << franja[h].Grp << endl;
	     h++ ;
	  }
	  if(h>maxovh) moco("la cagaste en la interface\n");
        } while( i != lirst[ibox] ); /*fin laso particulas del grid*/
      } /*fin iz*/
  } /*fin iy*/
 // Salida.close();
// cout << proc << " interface >> "<< ix<< endl ;
  return 0;
}

int change_coord(void)
{
  double cb ; 
  int i;
  cb = (boxsize/(double)(numprocs-1))*(double)proc - overhead ;

  for(i=0;i<particle.size();i++)
  {
    particle[i].Pos[0] = (float)((double)particle[i].Pos[0] - cb );
    if(particle[i].Pos[0] < 0.) 
	particle[i].Pos[0] =(float)((double)particle[i].Pos[0] + boxsize) ;
    if((double)particle[i].Pos[0] >= boxsize /*+2.*overhead*/) 
        particle[i].Pos[0] =(float)((double)particle[i].Pos[0] - boxsize) ;// - 2.*overhead; 

  }
  return 0;
}

// modificada para que sobre un micro
int read_header(void)
{
  int  npproc;
  int i,j,k;
  int numpartfile,numpartfilemax,numpart;
  char filename[200];
  FILE *fd;
  double boxsizet;
  
  printf("\033[22;33m proc: %d FILES: %d numprocs%d \033[22;0m \n",proc,FILES,numprocs);
  fflush(stdout);
  if( numprocs -1 < FILES )
  {
    BLUE; printf("proceso: %d\n",proc);
    moco("el numero de procesos es menor que el nunero de files\n");
    printf("numero de procesos: \033[22;32m %d \033[22;0m numero de files: \033[22;32m %d \033[22;0m \n",
            numprocs - 1 ,proc);  FIN;
  }

  numpartfile=0;  numpartfilemax=0;  numpart=0;
  if( proc < FILES )
  {
    sprintf(filename,"%dmpcdata/snapshot_0%d.%d",INTBOX,SNAPNUM,proc); 
//    cout<< filename << endl;
    fd=fopen(filename,"r"); if(fd==NULL) {cout << "@@@@"<<proc << endl;moco("fallo el open\n"); FIN;}

    skip(fd); my_fread(&header1, sizeof(header1), 1, fd); skip(fd); skip(fd);
  //  for(i=0;i<6;i++)printf("masa particula %f\n",header1.mass[i]);
  //  for(i=0;i<6;i++)printf("numero particula %d\n",header1.npart[i]);

    for(i=0;i<6;i++)
	    numpartfile+=header1.npart[i];
    if(proc==0)
	    for(i=0;i<6;i++){numpart+=header1.npartTotal[i];}
    
    boxsizet=(double)header1.BoxSize;
    if(proc==0)boxsize=(double)header1.BoxSize;
  }

  MPI_Bcast(&boxsize,1,MPI_DOUBLE,0,MPI_COMM_WORLD);

  if( proc < FILES && boxsizet!=boxsize)
  {
	  moco("Diferentes Boxsizes\n"); FIN;
  }

  MPI_Bcast(&numpart,1,MPI_INT,0,MPI_COMM_WORLD);
  nparticles_total=numpart;
  printf("\033[01;36m NUMERO TOTAL %d \033[01;0m\n",nparticles_total);fflush(stdout);

  numpartfilemax *= 1 + (int)ceil(2.*overhead/boxsize);
  if(proc < FILES )fclose(fd);
  return 0;
  
}

// modificada para que sobre un micro
int read_gadget_overlap(void)
{
  int npproc;
  int i,j,k,ss;
  int numpartfile,numpartfilemax,coco;
  char filename[200];
  FILE *fPos, *fVel, *fId;
  float xx[3],vv[3];
  int npaqs, ind;
  int nsend,nnn;
  float lolo;
  struct ptcl *posbuf;
  size_t buffersize;
  int contpos, npprocold;
  int ix,iy,iz;
  int iov,nov,id;
  double xdesp;
  double ladosize;
  int *npaq;
  struct ptcl *paq;
  struct ptcl tmptcl,nulaptcl;

  //  MPI_Datatype  Partmpi; es publico ahora
  MPI_Datatype  type[5]= {MPI_FLOAT, MPI_FLOAT, MPI_INT, MPI_INT, MPI::BOOL};
  int blocklen[5]      = { 3, 3, 1, 1, 1};
  MPI_Aint disp[5];
  
  numpartfile=0;  numpartfilemax=0;  

  //Abre el archivo lee el header e inicializa los punteros al file
  if( proc < FILES )
  {
    sprintf(filename,"%dmpcdata/snapshot_0%d.%d",INTBOX,SNAPNUM,proc); //printf(filename);

    fPos=fopen(filename,"r"); if(fPos==NULL){ moco("fallo el open\n"); FIN;}
    fVel=fopen(filename,"r"); if(fVel==NULL){ moco("fallo el open\n"); FIN;}
    fId =fopen(filename,"r"); if(fId ==NULL){ moco("fallo el open\n"); FIN;}

    fseek(fPos,sizeof(int),SEEK_CUR);
    my_fread(&header1, sizeof(header1), 1, fPos); 
//    for(i=0;i<6;i++)printf("masa particula %f\n",header1.mass[i]);
//    for(i=0;i<6;i++)printf("numero particula %d\n",header1.npart[i]);
    
    fseek(fPos,2*sizeof(int),SEEK_CUR);
    fseek(fVel,sizeof(header1)+3*sizeof(int),SEEK_CUR); 
    fseek(fId, sizeof(header1)+3*sizeof(int),SEEK_CUR); 

    for(i=0;i<6;i++)
 	    numpartfile+=header1.npart[i];

    fseek(fVel,sizeof(float)*3*numpartfile + 2*sizeof(int),SEEK_CUR); // para las vel saltamos un bloque
    fseek(fId, sizeof(float)*6*numpartfile + 4*sizeof(int),SEEK_CUR); // para los Ids 2 bloques
  }

  MPI_Allreduce(&numpartfile,&numpartfilemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD); 
//  BLUE;
//  printf("%d | %d \n",numpartfile,proc);
//  DEFA;

  ladosize = boxsize/((double)(numprocs-1));
  nov      = (int)floor(ladosize/overhead);
  overhead = ladosize/(double)nov; // redefine un overhead divisible
//  printf("overhead usado: %f boxisze %f\n",overhead,boxsize);

  numpartfilemax = numpartfilemax + (int)((double)numpartfilemax*(2.*overhead/boxsize));

  npaqs    = numpartfilemax/NPAQMAX +1; // numero de paquetes
  
  for(i=0;i<3;i++) nulaptcl.Pos[i] = 0. ;
  for(i=0;i<3;i++) nulaptcl.Vel[i] = 0. ;
                   nulaptcl.Id     = -1 ;
                   nulaptcl.Idgrupo   = -1 ;
                   nulaptcl.visitado  = false ;

 
  paq  =new struct ptcl [(numprocs-1)*NPAQMAX];
  npaq =new int [numprocs-1];


  // creamos la estructura para comunicar los paquetes en MPI
  disp[0] = (MPI_Aint)&paq[0].Pos    - (MPI_Aint)&paq[0];
  disp[1] = (MPI_Aint)&paq[0].Vel    - (MPI_Aint)&paq[0];
  disp[2] = (MPI_Aint)&paq[0].Id     - (MPI_Aint)&paq[0];
  disp[3] = (MPI_Aint)&paq[0].Idgrupo  - (MPI_Aint)&paq[0];
  disp[4] = (MPI_Aint)&paq[0].visitado  - (MPI_Aint)&paq[0];

  MPI_Type_create_struct(5, blocklen, disp, type, &Partmpi);
  MPI_Type_commit(&Partmpi);

  buffersize=(((NPAQMAX+MPI_BSEND_OVERHEAD)*sizeof(struct ptcl)));
  if(proc==0)
  {
    RED; printf("Buffer size is %f Mbytes\n",(float)buffersize/1024./1024.); 
    DEFA; fflush(stdout);
  }
  
  posbuf=(struct ptcl*)malloc((numprocs-1)*buffersize);

  contpos=0;
  MPI_Buffer_attach(posbuf,(numprocs-2)*buffersize);
  for(i=0;i<npaqs;i++)
  {
    for(j=0;j<(numprocs-1)*NPAQMAX;j++) 
	  paq[j]=nulaptcl; 
    if(proc==0)
    {
      printf(" %d Tranzando bloque: \033[01;36m %d / %d \033[01;0m\r",proc,i,npaqs);
      fflush(stdout);
    }

    for(j=0;j<numprocs-1;j++) npaq[j]=0;
    npproc=0;
      
    // lectura de los paquetes 
    if( proc < FILES )
    {
      for(j=0;j<NPAQMAX && contpos<numpartfile ;j++,contpos++)
      {
           my_fread(&xx,sizeof(float),3,fPos);
           my_fread(&vv,sizeof(float),3,fVel);
           my_fread(&id,sizeof(int  ),1,fId ); tmptcl.Id =id ;
	   //if(drand48()>0.5)continue;
           //controla periodicidad
	   tmptcl.Idgrupo=-1; // inicializandolo para la identificacion
	   tmptcl.visitado=false; // inicializandolo para la identificacion
           for( k = 0 ; k < 3 ; k++ )
           {
             if( xx[k] >= boxsize ) xx[k] -= boxsize;
             if( xx[k] < 0. )        xx[k] += boxsize;
             if( xx[k] < 0. ||  xx[k] >= boxsize  )  moco("-------b-size 3\n");
           }

       	   for(k=0;k<3;k++){ tmptcl.Pos[k]=xx[k] ;
                             tmptcl.Vel[k]=vv[k] ;}

           nsend=(int)floor((double)xx[0]/ladosize); //nodo al que va la ptcl

           if(nsend>= numprocs-1 || nsend < 0){moco("LA CAGASTE");printf("nsend: %d\n",nsend);FIN;}
           if(nsend==proc) npproc++; 
        
	   // descompone en rodajas disjuntas
	   paq[nsend*NPAQMAX + npaq[nsend]]=tmptcl; npaq[nsend]++;

           //aca agrega la superposicion 
           xdesp=ladosize*(double)nsend;
           iov=(int)floor((((double)xx[0]-xdesp)/ladosize)*(double)nov);
		if(iov <0 || iov>nov-1){ moco ("charola !!!!!!");FIN;}
           if(iov==0)
           {
               nsend--; if(nsend < 0) nsend=numprocs-2;
               if(nsend==proc) npproc++; 
	       paq[nsend*NPAQMAX + npaq[nsend]]=tmptcl; npaq[nsend]++;
           }
           if(iov==nov-1)
           {
               nsend++; if(nsend >= numprocs-1) nsend =0;
               if(nsend==proc) npproc++; 
	       paq[nsend*NPAQMAX + npaq[nsend]]=tmptcl; npaq[nsend]++;
           }
      } 

      //carga las particulas propias
      for(k=0;k<npproc;k++)
      {
         tmptcl = paq[proc*NPAQMAX + k];
	 particle.push_back(tmptcl) ;
      }

      for(j=0;j<numprocs-1;j++) 
      {
         if(j!=proc)
         {
           ind=j*NPAQMAX; 
	   MPI_Bsend(&paq[ind],NPAQMAX,Partmpi,j,33,MPI_COMM_WORLD);
	   for(k=0;k<npaq[j];k++) if(paq[ind+k].Id<0) {moco("SOS CULCPABLE");FIN;}
         }
      }
    } //fin del if de los lectores

    // lazo para recivir
    if( proc < numprocs - 1)
    {
	    for(j=0;j<FILES;j++)
	    {
		    if(j!=proc)
		    {
			    ind=j*NPAQMAX; 
		    	    MPI_Recv(&paq[ind],NPAQMAX,Partmpi,j,33,MPI_COMM_WORLD,&status);
			    try
			    {
			   	    for(k=0;k<NPAQMAX;k++)
			  		    if(paq[k+ind].Id>=0)
						    particle.push_back(paq[k+ind]) ;
			    }
		     	    catch(bad_alloc) 
			    {
			     	    cout << endl << particle.size()<< " ja " 
				       	    << sizeof(struct ptcl) << " " << proc <<flush << endl;
					    moco("bad_alloc!!");FIN;
			    }
		    }
	    }
    }

    if( proc < FILES )
    {
        k=0; 
	MPI_Buffer_detach(posbuf,&k);MPI_Buffer_attach(posbuf, k);
    }

  }

  //fin del lazo de los paquetes
  if( proc < FILES )
  {
    fclose(fPos);
    fclose(fVel);
    fclose(fId);
  }

//  printf("proc \033[22;31m%d \033[22;0m npproc %d numpartfilemax %d \n",proc,npproc,numpartfilemax);
  fflush(stdout);
  if(npproc>numpartfilemax)
  {
    moco("NPROC > NUMPARTFILEMAX: ");
    printf("proc=%d npproc=%d numpartfilemax=%d\n",proc,npproc,numpartfilemax);
    fflush(stdout);
    FIN;
  }

  free(posbuf);
  delete[] npaq; delete[] paq;
  return(particle.size());
/*Termino de leer y distribuir*/
}

int idenBFS(float r0)
{
	int i,j,u,v;
	float percent,pp;
	queue  <int> Q;
        vector <int> listvec;
	struct grp halotmp;

	cout << "I'm in!  " << proc << endl;
	fflush(stdout);
        if(proc==numprocs-1)return 1;
	ngrupos=0;
	pp=(float) particle.size();
	for(i=0;i<particle.size();i++)
	{   
            percent= ((float) i/pp)*100.;
	    if(fmodf(percent,10.)==0)
	    cout << "procesando "<< percent  << " % : proc: " << proc << endl;
            if(particle[i].visitado) continue ;
            particle[i].visitado  = true ;
            //Encolar(Q, s);
	    
	    halotmp.name = ngrupos;
	    halotmp.members.push_back(&particle[i]);
	    Q.push(i);
            do
            {
                  // extraemos el nodo u de la cola Q y exploramos 
	          // todos sus nodos adyacentes
                  u     = Q.front();
   	          Q.pop();
		  particle[u].Idgrupo= ngrupos;
		  busv(u,(float)r0,&listvec);
		  // devuelve la lista de vecinos no visitados
                  for(j=0;j<listvec.size();j++)
	          {
	               halotmp.members.push_back(&particle[listvec[j]]);
		       particle[listvec[j]].Idgrupo    = ngrupos;
                       if(particle[listvec[j]].visitado)
			       continue;
		       else
			       particle[listvec[j]].visitado = true;
		       Q.push(listvec[j]);  //Encolar v en Q
	          }
		  listvec.clear();
            }while(!Q.empty());//agotadas todas las 
	                       //adyacencias salimos
	    if(halotmp.members.size()>=2)
	    {
		    grupo.push_back(halotmp);
	            ngrupos++;
	    }
	    halotmp.members.clear();
	}
	cout << "I'm out!  " << proc << endl;
	return 0;
}

int grid(void)
{
  int i, j, k;
  int ix, iy, iz;

  //printf("grid\n");
  fflush(stdout);
  if(proc<numprocs -1)
  {
	  for( i = 0 ; i < nx ; i++)
	       	  for( j = 0 ; j < np ; j++)
		     	  for( k = 0 ; k < np ; k++)
			       	  lirst  [(i * np + j) * np + k] = -1 ;
  
	  for( i = 0 ; i < particle.size() ; i++ )  
		  ll[i] = -1 ;
	  //  for( i = 0 ; i < particle.size() ; i++ )
	  //  f(&particle[i].Pos[0]==NULL)moco("no data #*$?%!!!");
	  totalsizeL=0; totalsizeR=0;
	  //  cout  << proc  << "   "  << facbinx  << " ||| " <<  facbinp  << " *** "  << r0  << endl;
	  for( i = 0 ; i < particle.size() ; i++ )
	  {
	      	  ix=(int)floor(((double)particle[i].Pos[0])*facbinx)  ;
	      	  iy=(int)floor(((double)particle[i].Pos[1])*facbinp)  ;
	      	  iz=(int)floor(((double)particle[i].Pos[2])*facbinp)  ;
	      	  if(ix>=nx || ix < 0)
	      	  {
	      		  moco("---------------------la cagaste ix : ");
	      		  cout<< ix << " " 
		      		  << particle[i].Pos[0] << " " 
		      		  << facbinx << " " 
		      		  << facbinp << endl; 
	      	  }
	      	  if(iy>=np || iy < 0) moco("---------------------la cagaste iy");
	      	  if(iz>=np || iz < 0) moco("---------------------la cagaste iz");
	      	  lirst  [(ix * np + iy) * np + iz] = i ;
	      	  if(fabs(xsize-(double)particle[i].Pos[0]) < 2.*overhead )  totalsizeR++;
	      	  if((double)particle[i].Pos[0] < 2.*overhead )  totalsizeL++;
	  }  
  }else{ totalsizeR=0;totalsizeL=0;}

  cout << endl << proc << " proc, L "<<totalsizeL << " y R " <<  totalsizeR  << endl;
  fflush(stdout);
  MPI_Allreduce(&totalsizeR,&maxsizeR,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD); 
  MPI_Allreduce(&totalsizeL,&maxsizeL,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
  maxovh = maxsizeR > maxsizeL ? maxsizeR : maxsizeL;
//  printf("3-----\n");

  if(proc<numprocs-1)
  {
	  for( i = 0 ; i < particle.size() ; i++ )
	  {
		  ix=(int)floor((double)particle[i].Pos[0]*facbinx);
		  iy=(int)floor((double)particle[i].Pos[1]*facbinp);
	       	  iz=(int)floor((double)particle[i].Pos[2]*facbinp);
	       	  ll[lirst[(ix * np + iy) * np + iz]] = i ;
	      	  lirst[(ix * np + iy) * np + iz] = i ;
	  }  
  }

  printf("fin grid\n");
  fflush(stdout);
  return 0;
}

int busv(int ic, float r0, vector<int> *listvec)
{
  int ixc, iyc, izc;
  int ixci, iyci, izci;
  int ixcf, iycf, izcf;
  int ix, iy, iz;
  int ixx, iyy, izz;
  int ibox;
  float xx, yy, zz, dis;
  int i;

  //Esta busv es solo no periodica con grideado en rodajas

  ixc=(int)(particle[ic].Pos[0]/(float)xsize*nx);
  ixci=ixc-1;
  ixcf=ixc+1;
  iyc=(int)(particle[ic].Pos[1]/(float)boxsize*np);
  iyci=iyc-1;
  iycf=iyc+1;
  izc=(int)(particle[ic].Pos[2]/(float)boxsize*np);
  izci=izc-1;
  izcf=izc+1;

  if( ixci < 0 ) ixci = 0;
  if( iyci < 0 ) iyci = 0;
  if( izci < 0 ) izci = 0;
  if( ixcf >= nx ) ixcf = nx-1;
  if( iycf >= np ) iycf = np-1;
  if( izcf >= np ) izcf = np-1;

  for( ixx = ixci ; ixx <= ixcf ; ixx++)
  {
    ix=ixx;
    for( iyy = iyci ; iyy <= iycf ; iyy++)
    {
      iy=iyy;
      for( izz = izci ; izz <= izcf ; izz++)
      {
        iz=izz;
        ibox=(ix * np + iy) * np + iz ;
        i=lirst[ibox];
        if(i==-1)continue;
        do
        {
          i=ll[i];
	  if(particle[i].visitado && i != lirst[ibox] )continue;
	  if(particle[i].visitado && i == lirst[ibox] )break;
          xx=particle[i].Pos[0]-particle[ic].Pos[0];
          yy=particle[i].Pos[1]-particle[ic].Pos[1];
          zz=particle[i].Pos[2]-particle[ic].Pos[2];
          dis=sqrt(xx*xx+yy*yy+zz*zz);
          if( dis < (float)r0 && particle[i].Idgrupo == -1 && i != ic )
               (*listvec).push_back(i);

        } while( i != lirst[ibox] ); /*fin laso particulas del grid*/
      } /*fin izz*/
    } /*fin iyy*/
  } /*fin ixx*/
  return 0;
}

